function[u]=utility(rho,cons);
if rho==1
u=log(cons);
else
    u=1/(1-rho)*(cons.^(1-rho)-1);
end